39tn AIAA/ASMt/SAt/A5fcS Joint 
Propulsion Conference and Exhibit 

Huntsville, AL AIAA-2003-4919 

‘ 20-23 July 2003 


Calculation of Turbine Axial Thrust by Coupled CFD Simulations of the Main 
Flow Path and Secondary Cavity Flow in an SLI LOX Turbine 

D. J. Dorney 

NASA/George C. Marshall Space Flight Center 
Space Transportation Directorate 
MSFC, AL 35812 

Bogdan Marcu 
Ken Tran 
Scott Sargent 

The Boeing Company 
Rocketdyne Propulsion and Power 
6633 Canoga Avenue 
Canoga Park, CA 91309-7922 


ABSTRACT 

Each single reusable Space Launch Initiative 
(SLI) booster rocket is an engine operating at a 
record vacuum thrust level of over 730,000 Ibf 
using LOX and LH2. This thrust is more than 
10% greater than that of the Delta IV rocket, 
resulting in relatively large LOX and LH2 
turbopumps. Since the SLI rocket employs a 
staged combustion cycle the level of pressure is 
very high (thousands of psia). This high 
pressure creates many engineering challenges, 
including the balancing of axial, forces on the 
turbopumps. One of the main parameters in the 
calculation of the axial force is the cavity 
pressure upstream of the turbine disk. The flow 
in this cavity is very complex. The lack of 
understanding of this flow environment hinders 
the accurate prediction of axial thrust. In order to 
narrow down the uncertainty band around the 
actual turbine axial force, a coupled, unsteady 
computational methodology has been developed 
to simulate the interaction between the turbine 
main flow path and the cavity flow. The 
CORSAIR solver, an unsteady three- 
dimensional Navier-Stokes code for 
turbomachinery applications, was used to solve 
for both the main and the secondary flow fields. 
Turbine axial thrust values are presented in 
conjunction with the CFD simulation, together 
with several considerations regarding the turbine 
instrumentation for axial thrust estimations 
during test. 
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INTRODUCTION 

In the recent years the harnessing of increased 
computing power through parallelization of CFD 
codes has allowed turbomachinery engineers to 
extend their analyses to the secondary flow paths 
of turbomachines. The prediction of unsteady, 
three-dimensional viscous flow physics using CFD 
for the main compressor and turbine flow paths has 
been coupled with the flow through disk and 
housing cavities in an effort to better model and 
understand the actual flow in a turbomachine. 

One of the main motivations for the analysis of the 
secondary flow paths has been the necessity to 
understand and quantify the heat transfer and 
cooling capacities of such flows, as well as the 
performance impact the designer must pay in order 
to use such capacities. Early works in addressing 
cavity flows [1-4] was driven by the need to 
understand the cooling cavity flows in basic 
geometries. The early numerical results reported in 
the literature [such as Refs. 5 and 6, among an 
exhaustive list of valuable contributions] have been 
based on 2-D, steady flow models that do not 
completely account for the complex environment of 
rotor-stator interaction coupled with the cavity flow 
field. Recently [7], results of more complex 
simulations have been reported. Such simulations 
are based on the coupling of large and well- 
validated CFD solvers for simultaneous calculations 
of the turbine main flow path (e.g., NASA MS- 
TURBO) and the cavity flow field (e.g., NASA 
SCISEAL). Special attention is focused on grid 
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issues and solving the flow accurately though 
the seal packages. 

While a wealth of information pertaining to the 
flow of cooling streams in the disk cavity 
geometries typical to air-breathing gas turbine 
engines has been, and is being, delivered by the 
contributions of numerous authors, less data is 
available for the disk cavity geometries typically 
encountered in the turbopumps equipping liquid 
propellant rocket engines. Such geometries 
have the same level of complexity as the gas 
turbine engines, but a more severe operating 
environment. Typically, cryogenic turbopump 
turbines operate in a hot, high-pressure gas 
environment, separated by inches from the 
cryogenic fluid flow path of the pump. If the 
pumped fluid and the turbine working fluid are 
compatible, small quantities of cold vaporized 
fluid from the pump secondary flow paths may 
be allowed to leak into the turbine cavity. The 
cold net inflow mixes with hot gas in the cavity in 
a manner that is not currently understood. If the 
fluids are not compatible (for example, in a liquid 
oxygen turbopump powered by gaseous hot 
hydrogen) the turbine cavity and the pump 
secondary flows must be carefully organized to 
keep the fluids separated. One way to achieve 
separation is through an interstage seal that 
allows each leakage to drain out. The turbine 
cavity in this case has a net outflow. 

In the latter example of a design configuration, a 
significant uncertainty arises for the case of very 
large turbopumps; namely the magnitude of the 
axial thrust produced by the turbine disks. The 
energetic nature of pump operation requires a 
carefully organized balance between the axial 
loads produced on the pump and turbine rotating 
components. 

The present paper specifically addresses this 
uncertainty as related to the Space Launch 
Initiative RS-83 Engine Liquid Oxygen 
turbopump via coupled CFD simulations in 
which both the turbine main flow path and the 
cavity flow are simultaneously solved in the full 
three-dimensional, viscous unsteady regime. 
The cavity net flow, pressure distribution and 
resulting axial thrust are calculated for three 
situations, corresponding to a well sealed cavity 
with minimum mass flow leaking into the drain 
system, and two cavities with increased levels of 
wear for the seals. 


LOX TURBINE AERODYNAMIC DESIGN 

The oxidizer turbopump (OTP) turbine for the 
Boeing Rocketdyne RS-83 hydrogen fueled engine 
is a single stage, subsonic design. The primary 
design objectives for this turbine are high 
performance and robust design. To satisfy these 
conditions the turbine achieves an efficiency of 
77% at the design pressure ratio (PR) of 1.22. 
Turbine design total inlet pressure is 3876 psia. 
The turbine has a relatively large pitch line diameter 
of 16.21 in. with a blade height to diameter ratio 
(L/D) of 8.9% and a throat aspect ratio (L/O) of 
4.69. While the pitch diameter of the turbine is 
large, the operation of the seal package preventing 
excessive mass flow leakage out of the turbine 
cavity imposes a low diameter size for the seal. 
Consequently, the turbine disk offers a large area 
for the cavity pressure to act upon and produce 
axial load. The understanding of the cavity flow and 
the associated pressure distribution becomes a 
strong requirement for a robust design. 

NUMERICAL PROCEDURE 

The governing equations considered in this study 
are the time-dependent, three-dimensional 
Reynolds-averaged Navier-Stokes equations. The 
numerical algorithm used in the computational 
procedure consists of a time marching, implicit, 
finite-difference scheme. The procedure is third- 
order spatially accurate and second-order 
temporally accurate. The inviscid fluxes are 
discretized according to the scheme developed by 
Roe [8]. The viscous fluxes are calculated using 
standard central differences. An approximate- 
factorization technique is used to compute the time 
rate changes in the primary variables. Newton sub- 
iterations are used at each global time step to 
increase stability and reduce linearization errors. 
The equations of motion are extended to turbulent 
flows using an eddy viscosity formulation. The 
turbulent viscosity is calculated using the two-layer 
Baldwin-Lomax algebraic turbulence model [9]. 

Message Passing Interface (MPI) and OpenMP 
directives have been implemented into the code to 
reduce the computation time for large-scale three- 
dimensional simulations. The use of MPI allows the 
coupling of different geometric components, such 
as the turbine main flow path and cavity, in a 
straightforward manner. This procedure was 
previously demonstrated for unsteady flow in a 
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supersonic turbine with straight centerline 
nozzles [10]. 

The computational procedure uses O- and H- 
type zonal grids to discretize the flow field and 
facilitate relative motion of the rotating 
components. The O-grids are body-fitted to the 
surfaces of the airfoils and generated using an 
elliptic equation solution procedure. They are 
used to properly resolve the viscous flow in the 
blade passages and to easily apply the algebraic 
turbulence model. The algebraically-generated 
H-grids are used to discretize the remainder of 
the flow field, including the turbine cavity. 

Figure 1 shows an axial cut through the 
computational grids. The turbine is modeled as a 
53 nozzle/53 blade turbine (1/1 computationally) 
using 111x45x31 (H) and 121x45x31 (O) grids 
for nozzle, and 111x45x31 (H) and 121x45x31 
(O) grids for rotor. The cavity is modeled using a 
35x45x111 (H) grid. 

The circumferential span of the computational 
domain is 6.8 °. 

The computational analysis has been validated 
on several rocket-engine turbine geometries, 
including coupled components such as straight 
centerline nozzles (e.g., Refs. 10-14). 

Boundary Conditions 

The theory of characteristics is used to 
determine the boundary conditions at the inlet 
and exit of the computational domain. For 
subsonic inlet flow the total pressure, total 
temperature, and the circumferential and radial 
flow angles are specified as a function of the 
radius. The upstream running Riemann 
invariant, is extrapolated from the interior of the 
computational domain. 

For subsonic outflow the circumferential and 
radial flow angles, total pressure, and the total 
temperature are extrapolated from the interior of 
the computational domain. For the turbine main 
flow path the total-to-static pressure ratio is 
specified at mid-span of the computational exit 
and the pressure at all other radial locations at 
the exit is obtained by integrating the equation 
for radial equilibrium. For the cavity flow, the 
upper computational domain is coupled with the 
main flow computational domain at the H-grid 
area corresponding to the turbine nozzle hub 
wall. The cavity walls are one stationary and one 


rotating with so slip boundary conditions enforced 
on each. 

Given the relatively large size of the grid, the 
labyrinth seal is not modeled numerically, as such 
additional modeling would have rendered the study 
computationally too expensive. Instead, the degree 
of wear in the seal is modeled through a pressure 
drop across the discharge area at the bottom of the 
cavity, with the largest the pressure drop 
corresponding to the most extensive wear modeled. 
The objective of the calculation is to model the 
cavity flow at progressively larger mass flow rates. 

Periodicity is enforced along the outer boundaries 
of the H-grids in the circumferential direction. For 
viscous simulations, no-slip boundary conditions 
are enforced along the solid surfaces. It is 
assumed that the normal derivative of the pressure 
is zero at solid wall surfaces. In addition, a 
specified (zero) heat flux is held constant in time 
along the solid surfaces. 

Figure 2 shows an instantaneous snapshot of the 
confluence of the two computational domains and a 
vector flow visualization of the flow pattern as mass 
flow is captured by the cavity volume. 

RESULTS 

The results presented in this report are 
circumferentially and time averaged unless 
specified otherwise. Of interest is the radial 
distribution of pressures p=p(r) for the three 
different flow regimes considered. Table 1 
presents the mass flow values for both the cavity 
and the turbine flow path for the three regimes 
denoted as Low Flow (which is the nominal flow as 
per design), Medium Flow and High Flow. All 
denoted values refer to the mass flow leaked at the 
bottom of the cavity. The algebraic sum of mass 
flows is not conserved from one case to another: 
the turbine boundary conditions, which are inlet 
pressure and temperature and inlet total to exit 
static pressure ratio, together with the coupling of 
the secondary cavity flow generate slightly different 
flow mass values entering the turbine nozzle. Thus, 
there is some overall variation in the turbine mass 
flow between the different cases. The variation, 
however, is less than one tenth of one percent. 
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Table 1 Mass flow values for the three cases 
(Ibm/s) 


Case 

Low Flow 

Med Flow 

High Flow 

Turbine 

257.4 

256.7 

256.9 

Cavity 

3.5 

4.8 

8.1 


section is modeled here) in relation to the rotor- 
stator instantaneous position may reduce the 
percentage variations shown in Table 2. However, 
no such correlation has been identified, the 
pressure variations shown in Fig. 5 being a 
combined effect of variations in pressure at the 
upper lip of the cavity and compressibility effects. 


Figure 3 shows the radial pressure distribution in 
the cavity for the three cavity mass flows. The 
average cavity pressure decreases proportional 
to the pressure drop across the discharge area 
at the bottom of the cavity and increasing mass 
flow. The shape of the radial pressure 
distribution curve also changes slightly. The fast 
pressure drop at the top of the cavity for all three 
cases is explained by examining the vector plot 
of the cavity flow field in Fig. 4. At the upper lip 
of the cavity the entering flow forms a stagnation 
area followed by a downward accelerated flow 
along the rotating disk wall. The “hump” in all 
three curves at a radius between 5.75 inches 
and 6 inches is due to the local geometry (locally 
increased volume) of the cavity. 

Table 2 lists the axial loads produced on the 
turbine disk by the three pressure distributions 
shown in Fig. 3. The values are calculated by 

R max 

F = 2k \p{r)rdr - (1 ) 

/firdn 

where p(r) is the radial pressure distribution, R min 
and Rmax are the radii at the top and bottom of 
the cavity, and P exit is the turbine discharge 
pressure 


Table 2. Axial load on the turbine disk 
(values in Ibf) 


Case 

Low Flow 

Med Flow 

High Flow 

F 

-17904.1 

-19121.7 

-21751.9 

F % 

variation 

±3.2% 

±3.4% 

±4.9% 


The percent variation in axial thrust is due to 
significant fluctuation in the pressure distribution 
values for each case. Figures 5a, 5b and 5c 
show each of the p = p(r) curves from Fig. 3, as 
well as it’s fluctuation range. The authors did not 
identify a strong correlation between the rotor 
blade position and the pressure fluctuation. If 
such a correlation existed, clocking effects for 
each cavity sector (of which only a 6.8° period 


The pressure distribution p = p(r) is never available 
when testing hardware, the price and complexity of 
the instrumentation would be prohibitive. A 
common practice is to model the pressure variation 
in the cavity based on a swirling flow model and 
compute a pressure distribution anchored on the 
pressure value at the lip of the cavity pfr^, which 
is known at design time from the turbine mean flow 
path analysis: 


P( r ™x)-P(r) 


k 2 p ^ 2 (rl~r 2 ) 
2gl44 


( 2 ) 


where p is the fluid density, £2 is the turbine angular 
velocity (in radians per second), g is the 
gravitational constant, and k is what is usually 
called as the pumping coefficient or slip factor 
coefficient 



The pumping coefficient is a way to estimate how 
much of the disk tangential velocity U is imparted to 
the fluid velocity, Cu. 

The true values of the k factor at each radius k=k(r) 
are shown in Fig. 6. If one were to use these 
values (smoothed) and compute the radial pressure 
distribution in the cavity using expression (2), the 
results would look as shown in Fig. 7. Clearly, at 
least for the geometry analyzed here, the swirling 
flow model is not quite accurate. 

In practice, for turbopump development testing one 
or several redundant pressure sensors are placed 
in the cavity at one particular radius. The choice of 
where to place the sensors depends on the practice 
of the team, the availability of space to place the 
sensor itself, or other particular issues. Figure 8 
shows the cavity analyzed in this study and four 
radial locations at which a static pressure sensor 
may be placed. These are four distinct possibilities. 
For each of these test configurations, the local 
pressure p(rprobe) value is used to calculate a value 
for the k factor using Eqn. (2). Then, also using 
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Eqn. (2) and the freshly calculated k factor, a 
radial pressure distribution is calculated, and 
integrated to obtain the axial disk load. Figure 9 
shows the comparison between the radial 
pressure distributions shown in Fig. 3 and the 
pressure distributions obtained based on the 
four static pressure measurement points, for 
each flow case. With the exception of the top 
portion, the radial pressure distribution is 
matched well for the cases when the choice is to 
place the pressure sensor at mid-to-lower radii 
rather than higher radii in the cavity. Optimum 
for the present case is the KP2 location at the 
5.75 in radius, which is also the location of the 
bulkier local volume in the cavity geometry. This 
observation is valid for the cavity geometry at 
hand, but a generalization cannot be made 
without analyzing a large number of other 
geometries of cavities. Table 3 lists the values 
for the axial turbine disk loads computed for the 
pressure distributions shown in Fig. 9 using 
expression (1). 


Table 3. Turbine disk axial loads (Ibf). 


Case 

Low Flow 

Med Flow 

High Flow 

F(lbf) 

%err 

F(lbf) 

%err 

F(lbf) 

%err 

KP1 

-19785 

10.5 

-20957 

9.6 

-24538 

12.8 

KP2 

-17611 

1.6 

-19106 

0.1 

-21282 

2.1 

KP3 

-17485 

2.3 

-18467 

3.4 

-20660 

5.0 

KP4 

-17185 

4.0 

-18260 

4.5 - 

-20495 

5.8 
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A series of numerical simulations have been 
performed for a coupled turbine/cavity geometry. 
The simulations were performed for three 
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predicted results differed from traditional swirl 
models, but provide valuable insight into the 
production of the turbine disk axial thrust loads. 
The application of a CFD analysis to determine 
sensor placement locations in turbine cavities 
was demonstrated. 
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Figure 2. Interface between the turbine main flow path and cavity flow: fluid is being ingested into the 

cavity. 
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Figure 3. Radial pressure distribution in the cavity. 
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Figure 5. Cavity radial pressure distributions and the associate ranges of variations for a) Low Flow, b) 
Medium Flow and c) High Flow. 
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Figure 6. Radial distribution of K factor values. K = Cu/U. 
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Figure 7. Comparison between cavity radial pressure distributions as computed from the CFD flow 
solution and from the swirl flow model. The K factor values from figure 6 have been used, after 

smoothing. 
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Figure 8. Pressure sensor probe locations. 
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Figure 9. Cavity radial pressure distribution calculated based on the various K factors anchored based on 

pressure measurement at the four KP points shown in figure 8. 
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